﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace LowLevelGraphics.Filter.Segmentations
{
    // SLIC.cpp: implementation of the SLIC class.
    //
    // Copyright (C) Radhakrishna Achanta 2012
    // All rights reserved
    // Email: firstname.lastname@epfl.ch
    //////////////////////////////////////////////////////////////////////

    //////////////////////////////////////////////////////////////////////
    // Construction/Destruction
    //////////////////////////////////////////////////////////////////////

    unsafe public class SLIC
    {
        int m_width;
        int m_height;
        int m_depth;

        double* m_lvec;
        double* m_avec;
        double* m_bvec;

        double** m_lvecvec;
        double** m_avecvec;
        double** m_bvecvec;

        public SLIC()
        {
            m_lvec = null;
            m_avec = null;
            m_bvec = null;

            m_lvecvec = null;
            m_avecvec = null;
            m_bvecvec = null;
        }

        ~SLIC()
        {
            //if(m_lvec) delete [] m_lvec;
            //if(m_avec) delete [] m_avec;
            //if(m_bvec) delete [] m_bvec;


            //if(m_lvecvec)
            //{
            //    for( int d = 0; d < m_depth; d++ ) delete [] m_lvecvec[d];
            //    delete [] m_lvecvec;
            //}
            //if(m_avecvec)
            //{
            //    for( int d = 0; d < m_depth; d++ ) delete [] m_avecvec[d];
            //    delete [] m_avecvec;
            //}
            //if(m_bvecvec)
            //{
            //    for( int d = 0; d < m_depth; d++ ) delete [] m_bvecvec[d];
            //    delete [] m_bvecvec;
            //}
        }

        //==============================================================================
        ///	RGB2XYZ
        ///
        /// sRGB (D65 illuninant assumption) to XYZ conversion
        //==============================================================================
        public void RGB2XYZ(
            int sR,
            int sG,
            int sB,
            double X,
            double Y,
            double Z)
        {
            double R = sR / 255.0;
            double G = sG / 255.0;
            double B = sB / 255.0;

            double r, g, b;

            if (R <= 0.04045) r = R / 12.92;
            else r = Math.Pow((R + 0.055) / 1.055, 2.4);
            if (G <= 0.04045) g = G / 12.92;
            else g = Math.Pow((G + 0.055) / 1.055, 2.4);
            if (B <= 0.04045) b = B / 12.92;
            else b = Math.Pow((B + 0.055) / 1.055, 2.4);

            X = r * 0.4124564 + g * 0.3575761 + b * 0.1804375;
            Y = r * 0.2126729 + g * 0.7151522 + b * 0.0721750;
            Z = r * 0.0193339 + g * 0.1191920 + b * 0.9503041;
        }

        //===========================================================================
        ///	RGB2LAB
        //===========================================================================
        public void RGB2LAB(int sR, int sG, int sB, double lval, double aval, double bval)
        {
            //------------------------
            // sRGB to XYZ conversion
            //------------------------
            double X = 0, Y = 0, Z = 0;
            RGB2XYZ(sR, sG, sB, X, Y, Z);

            //------------------------
            // XYZ to LAB conversion
            //------------------------
            double epsilon = 0.008856;	//actual CIE standard
            double kappa = 903.3;		//actual CIE standard

            double Xr = 0.950456;	//reference white
            double Yr = 1.0;		//reference white
            double Zr = 1.088754;	//reference white

            double xr = X / Xr;
            double yr = Y / Yr;
            double zr = Z / Zr;

            double fx, fy, fz;
            if (xr > epsilon) fx = Math.Pow(xr, 1.0 / 3.0);
            else fx = (kappa * xr + 16.0) / 116.0;
            if (yr > epsilon) fy = Math.Pow(yr, 1.0 / 3.0);
            else fy = (kappa * yr + 16.0) / 116.0;
            if (zr > epsilon) fz = Math.Pow(zr, 1.0 / 3.0);
            else fz = (kappa * zr + 16.0) / 116.0;

            lval = 116.0 * fy - 16.0;
            aval = 500.0 * (fx - fy);
            bval = 200.0 * (fy - fz);
        }

        //===========================================================================
        ///	DoRGBtoLABConversion
        ///
        ///	For whole image: overlaoded floating point version
        //===========================================================================
        public void DoRGBtoLABConversion(
            uint* ubuff,
            double[] lvec,
            double[] avec,
            double[] bvec)
        {
            int sz = m_width * m_height;
            lvec = new double[sz];
            avec = new double[sz];
            bvec = new double[sz];

            for (int j = 0; j < sz; j++)
            {
                int r = (int)(ubuff[j] >> 16) & 0xFF;
                int g = (int)(ubuff[j] >> 8) & 0xFF;
                int b = (int)(ubuff[j]) & 0xFF;

                RGB2LAB(r, g, b, lvec[j], avec[j], bvec[j]);
            }
        }

        //===========================================================================
        ///	DoRGBtoLABConversion
        ///
        /// For whole volume
        //===========================================================================
        public void DoRGBtoLABConversion(
            uint** ubuff,
            double** lvec,
            double** avec,
            double** bvec)
        {
            int sz = m_width * m_height;
            for (int d = 0; d < m_depth; d++)
            {
                for (int j = 0; j < sz; j++)
                {
                    int r = (int)(ubuff[d][j] >> 16) & 0xFF;
                    int g = (int)(ubuff[d][j] >> 8) & 0xFF;
                    int b = (int)(ubuff[d][j]) & 0xFF;

                    RGB2LAB(r, g, b, lvec[d][j], avec[d][j], bvec[d][j]);
                }
            }
        }

        //=================================================================================
        /// DrawContoursAroundSegments
        ///
        /// Internal contour drawing option exists. One only needs to comment the if
        /// statement inside the loop that looks at neighbourhood.
        //=================================================================================
        public void DrawContoursAroundSegments(
            uint* ubuff,
            int* labels,
            int width,
            int height,
            uint color)
        {
            int[] dx8 = new int[] { -1, -1, 0, 1, 1, 1, 0, -1 };
            int[] dy8 = new int[] { 0, -1, -1, -1, 0, 1, 1, 1 };

            	int sz = width*height;

            /*    vector<bool> istaken(sz, false);

                int mainindex(0);
                for( int j = 0; j < height; j++ )
                {
                    for( int k = 0; k < width; k++ )
                    {
                        int np(0);
                        for( int i = 0; i < 8; i++ )
                        {
                            int x = k + dx8[i];
                            int y = j + dy8[i];

                            if( (x >= 0 && x < width) && (y >= 0 && y < height) )
                            {
                                int index = y*width + x;

                                if( false == istaken[index] )//comment this to obtain internal contours
                                {
                                    if( labels[mainindex] != labels[index] ) np++;
                                }
                            }
                        }
                        if( np > 1 )//change to 2 or 3 for thinner lines
                        {
                            ubuff[mainindex] = color;
                            istaken[mainindex] = true;
                        }
                        mainindex++;
                    }
                }*/


            //int sz = width*height;
            //vector<bool> istaken(sz, false);
            List<int> contourx = new List<int>(sz); List<int> contoury = new List<int>(sz);
            List<int> mainindex = new List<int>(0);List<int> cind = new List<int>(0);
            //for( int j = 0; j < height; j++ )
            //{
            //    for( int k = 0; k < width; k++ )
            //    {
            //        int np(0);
            //        for( int i = 0; i < 8; i++ )
            //        {
            //            int x = k + dx8[i];
            //            int y = j + dy8[i];

            //            if( (x >= 0 && x < width) && (y >= 0 && y < height) )
            //            {
            //                int index = y*width + x;

            //                //if( false == istaken[index] )//comment this to obtain internal contours
            //                {
            //                    if( labels[mainindex] != labels[index] ) np++;
            //                }
            //            }
            //        }
            //        if( np > 1 )
            //        {
            //            contourx[cind] = k;
            //            contoury[cind] = j;
            //            istaken[mainindex] = true;
            //            //img[mainindex] = color;
            //            cind++;
            //        }
            //        mainindex++;
            //    }
            //}

            int numboundpix = cind.Count;// int numboundpix = cind//int(contourx.size());
            for (int j = 0; j < numboundpix; j++)
            {
                int ii = contoury[j] * width + contourx[j];
                ubuff[ii] = 0xffffff;

                for (int n = 0; n < 8; n++)
                {
                    int x = contourx[j] + dx8[n];
                    int y = contoury[j] + dy8[n];
                    if ((x >= 0 && x < width) && (y >= 0 && y < height))
                    {
                        int ind = y * width + x;
                        //if(!istaken[ind]) ubuff[ind] = 0;
                    }
                }
            }
        }


        //==============================================================================
        ///	DetectLabEdges
        //==============================================================================
        //public void DetectLabEdges(
        //    double*				lvec,
        //    double*				avec,
        //    double*				bvec,
        //    int					width,
        //    int					height,
        //    vector<double>&				edges)
        //{
        //    int sz = width*height;

        //    //edges.resize(sz,0);
        //    //for( int j = 1; j < height-1; j++ )
        //    //{
        //    //    for( int k = 1; k < width-1; k++ )
        //    //    {
        //    //        int i = j*width+k;

        //    //        double dx = (lvec[i-1]-lvec[i+1])*(lvec[i-1]-lvec[i+1]) +
        //    //                    (avec[i-1]-avec[i+1])*(avec[i-1]-avec[i+1]) +
        //    //                    (bvec[i-1]-bvec[i+1])*(bvec[i-1]-bvec[i+1]);

        //    //        double dy = (lvec[i-width]-lvec[i+width])*(lvec[i-width]-lvec[i+width]) +
        //    //                    (avec[i-width]-avec[i+width])*(avec[i-width]-avec[i+width]) +
        //    //                    (bvec[i-width]-bvec[i+width])*(bvec[i-width]-bvec[i+width]);

        //    //        //edges[i] = fabs(dx) + fabs(dy);
        //    //        edges[i] = dx*dx + dy*dy;
        //    //    }
        //    //}
        //}

        //===========================================================================
        ///	PerturbSeeds
        //===========================================================================
        //public void PerturbSeeds(
        //    vector<double>&				kseedsl,
        //    vector<double>&				kseedsa,
        //    vector<double>&				kseedsb,
        //    vector<double>&				kseedsx,
        //    vector<double>&				kseedsy,
        //        const vector<double>&                   edges)
        //{
        //    int[] dx8 = new int[]{-1, -1,  0,  1, 1, 1, 0, -1};
        //    int[] dy8 = new int[]{ 0, -1, -1, -1, 0, 1, 1,  1};

        //    int numseeds = kseedsl.size();

        //    for( int n = 0; n < numseeds; n++ )
        //    {
        //        int ox = kseedsx[n];//original x
        //        int oy = kseedsy[n];//original y
        //        int oind = oy*m_width + ox;

        //        int storeind = oind;
        //        for( int i = 0; i < 8; i++ )
        //        {
        //            int nx = ox+dx8[i];//new x
        //            int ny = oy+dy8[i];//new y

        //            if( nx >= 0 && nx < m_width && ny >= 0 && ny < m_height)
        //            {
        //                int nind = ny*m_width + nx;
        //                if( edges[nind] < edges[storeind])
        //                {
        //                    storeind = nind;
        //                }
        //            }
        //        }
        //        if(storeind != oind)
        //        {
        //            kseedsx[n] = storeind%m_width;
        //            kseedsy[n] = storeind/m_width;
        //            kseedsl[n] = m_lvec[storeind];
        //            kseedsa[n] = m_avec[storeind];
        //            kseedsb[n] = m_bvec[storeind];
        //        }
        //    }
        //}


        //===========================================================================
        ///	GetLABXYSeeds_ForGivenStepSize
        ///
        /// The k seed values are taken as uniform spatial pixel samples.
        //===========================================================================
        //public void GetLABXYSeeds_ForGivenStepSize(
        //    vector<double>&				kseedsl,
        //    vector<double>&				kseedsa,
        //    vector<double>&				kseedsb,
        //    vector<double>&				kseedsx,
        //    vector<double>&				kseedsy,
        //    const int&					STEP,
        //    const bool&					perturbseeds,
        //    const vector<double>&       edgemag)
        //{
        //    const bool hexgrid = false;
        //    int numseeds(0);
        //    int n(0);

        //    //int xstrips = m_width/STEP;
        //    //int ystrips = m_height/STEP;
        //    int xstrips = (0.5+double(m_width)/double(STEP));
        //    int ystrips = (0.5+double(m_height)/double(STEP));

        //    int xerr = m_width  - STEP*xstrips;if(xerr < 0){xstrips--;xerr = m_width - STEP*xstrips;}
        //    int yerr = m_height - STEP*ystrips;if(yerr < 0){ystrips--;yerr = m_height- STEP*ystrips;}

        //    double xerrperstrip = double(xerr)/double(xstrips);
        //    double yerrperstrip = double(yerr)/double(ystrips);

        //    int xoff = STEP/2;
        //    int yoff = STEP/2;
        //    //-------------------------
        //    numseeds = xstrips*ystrips;
        //    //-------------------------
        //    kseedsl.resize(numseeds);
        //    kseedsa.resize(numseeds);
        //    kseedsb.resize(numseeds);
        //    kseedsx.resize(numseeds);
        //    kseedsy.resize(numseeds);

        //    for( int y = 0; y < ystrips; y++ )
        //    {
        //        int ye = y*yerrperstrip;
        //        for( int x = 0; x < xstrips; x++ )
        //        {
        //            int xe = x*xerrperstrip;
        //            int seedx = (x*STEP+xoff+xe);
        //            if(hexgrid){ seedx = x*STEP+(xoff<<(y&0x1))+xe; seedx = min(m_width-1,seedx); }//for hex grid sampling
        //            int seedy = (y*STEP+yoff+ye);
        //            int i = seedy*m_width + seedx;

        //            kseedsl[n] = m_lvec[i];
        //            kseedsa[n] = m_avec[i];
        //            kseedsb[n] = m_bvec[i];
        //            kseedsx[n] = seedx;
        //            kseedsy[n] = seedy;
        //            n++;
        //        }
        //    }


        //    if(perturbseeds)
        //    {
        //        PerturbSeeds(kseedsl, kseedsa, kseedsb, kseedsx, kseedsy, edgemag);
        //    }
        //}

        //===========================================================================
        ///	GetKValues_LABXYZ
        ///
        /// The k seed values are taken as uniform spatial pixel samples.
        //===========================================================================
        //public void GetKValues_LABXYZ(
        //    vector<double>&				kseedsl,
        //    vector<double>&				kseedsa,
        //    vector<double>&				kseedsb,
        //    vector<double>&				kseedsx,
        //    vector<double>&				kseedsy,
        //    vector<double>&				kseedsz,
        //        const int&				STEP)
        //{
        //    const bool hexgrid = false;
        //    int numseeds(0);
        //    int n(0);

        //    int xstrips = (0.5+double(m_width)/double(STEP));
        //    int ystrips = (0.5+double(m_height)/double(STEP));
        //    int zstrips = (0.5+double(m_depth)/double(STEP));

        //    int xerr = m_width  - STEP*xstrips;if(xerr < 0){xstrips--;xerr = m_width - STEP*xstrips;}
        //    int yerr = m_height - STEP*ystrips;if(yerr < 0){ystrips--;yerr = m_height- STEP*ystrips;}
        //    int zerr = m_depth  - STEP*zstrips;if(zerr < 0){zstrips--;zerr = m_depth - STEP*zstrips;}

        //    double xerrperstrip = double(xerr)/double(xstrips);
        //    double yerrperstrip = double(yerr)/double(ystrips);
        //    double zerrperstrip = double(zerr)/double(zstrips);

        //    int xoff = STEP/2;
        //    int yoff = STEP/2;
        //    int zoff = STEP/2;
        //    //-------------------------
        //    numseeds = xstrips*ystrips*zstrips;
        //    //-------------------------
        //    kseedsl.resize(numseeds);
        //    kseedsa.resize(numseeds);
        //    kseedsb.resize(numseeds);
        //    kseedsx.resize(numseeds);
        //    kseedsy.resize(numseeds);
        //    kseedsz.resize(numseeds);

        //    for( int z = 0; z < zstrips; z++ )
        //    {
        //        int ze = z*zerrperstrip;
        //        int d = (z*STEP+zoff+ze);
        //        for( int y = 0; y < ystrips; y++ )
        //        {
        //            int ye = y*yerrperstrip;
        //            for( int x = 0; x < xstrips; x++ )
        //            {
        //                int xe = x*xerrperstrip;
        //                int i = (y*STEP+yoff+ye)*m_width + (x*STEP+xoff+xe);

        //                kseedsl[n] = m_lvecvec[d][i];
        //                kseedsa[n] = m_avecvec[d][i];
        //                kseedsb[n] = m_bvecvec[d][i];
        //                kseedsx[n] = (x*STEP+xoff+xe);
        //                kseedsy[n] = (y*STEP+yoff+ye);
        //                kseedsz[n] = d;
        //                n++;
        //            }
        //        }
        //    }
        //}

        //===========================================================================
        ///	PerformSuperpixelSLIC
        ///
        ///	Performs k mean segmentation. It is fast because it looks locally, not
        /// over the entire image.
        //===========================================================================
        //public void PerformSuperpixelSLIC(
        //    vector<double>&				kseedsl,
        //    vector<double>&				kseedsa,
        //    vector<double>&				kseedsb,
        //    vector<double>&				kseedsx,
        //    vector<double>&				kseedsy,
        //        int*&					klabels,
        //        const int&				STEP,
        //        const vector<double>&                   edgemag,
        //    const double&				M)
        //{
        //    int sz = m_width*m_height;
        //    const int numk = kseedsl.size();
        //    //----------------
        //    int offset = STEP;
        //        //if(STEP < 8) offset = STEP*1.5;//to prevent a crash due to a very small step size
        //    //----------------

        //    vector<double> clustersize(numk, 0);
        //    vector<double> inv(numk, 0);//to store 1/clustersize[k] values

        //    vector<double> sigmal(numk, 0);
        //    vector<double> sigmaa(numk, 0);
        //    vector<double> sigmab(numk, 0);
        //    vector<double> sigmax(numk, 0);
        //    vector<double> sigmay(numk, 0);
        //    vector<double> distvec(sz, DBL_MAX);

        //    double invwt = 1.0/((STEP/M)*(STEP/M));

        //    int x1, y1, x2, y2;
        //    double l, a, b;
        //    double dist;
        //    double distxy;
        //    for( int itr = 0; itr < 10; itr++ )
        //    {
        //        distvec.assign(sz, DBL_MAX);
        //        for( int n = 0; n < numk; n++ )
        //        {
        //                        y1 = max(0.0,			kseedsy[n]-offset);
        //                        y2 = min((double)m_height,	kseedsy[n]+offset);
        //                        x1 = max(0.0,			kseedsx[n]-offset);
        //                        x2 = min((double)m_width,	kseedsx[n]+offset);


        //            for( int y = y1; y < y2; y++ )
        //            {
        //                for( int x = x1; x < x2; x++ )
        //                {
        //                    int i = y*m_width + x;

        //                    l = m_lvec[i];
        //                    a = m_avec[i];
        //                    b = m_bvec[i];

        //                    dist =			(l - kseedsl[n])*(l - kseedsl[n]) +
        //                                    (a - kseedsa[n])*(a - kseedsa[n]) +
        //                                    (b - kseedsb[n])*(b - kseedsb[n]);

        //                    distxy =		(x - kseedsx[n])*(x - kseedsx[n]) +
        //                                    (y - kseedsy[n])*(y - kseedsy[n]);

        //                    //------------------------------------------------------------------------
        //                    dist += distxy*invwt;//dist = sqrt(dist) + sqrt(distxy*invwt);//this is more exact
        //                    //------------------------------------------------------------------------
        //                    if( dist < distvec[i] )
        //                    {
        //                        distvec[i] = dist;
        //                        klabels[i]  = n;
        //                    }
        //                }
        //            }
        //        }
        //        //-----------------------------------------------------------------
        //        // Recalculate the centroid and store in the seed values
        //        //-----------------------------------------------------------------
        //        //instead of reassigning memory on each iteration, just reset.

        //        sigmal.assign(numk, 0);
        //        sigmaa.assign(numk, 0);
        //        sigmab.assign(numk, 0);
        //        sigmax.assign(numk, 0);
        //        sigmay.assign(numk, 0);
        //        clustersize.assign(numk, 0);
        //        //------------------------------------
        //        //edgesum.assign(numk, 0);
        //        //------------------------------------

        //        {int ind(0);
        //        for( int r = 0; r < m_height; r++ )
        //        {
        //            for( int c = 0; c < m_width; c++ )
        //            {
        //                sigmal[klabels[ind]] += m_lvec[ind];
        //                sigmaa[klabels[ind]] += m_avec[ind];
        //                sigmab[klabels[ind]] += m_bvec[ind];
        //                sigmax[klabels[ind]] += c;
        //                sigmay[klabels[ind]] += r;
        //                //------------------------------------
        //                //edgesum[klabels[ind]] += edgemag[ind];
        //                //------------------------------------
        //                clustersize[klabels[ind]] += 1.0;
        //                ind++;
        //            }
        //        }}

        //        {for( int k = 0; k < numk; k++ )
        //        {
        //            if( clustersize[k] <= 0 ) clustersize[k] = 1;
        //            inv[k] = 1.0/clustersize[k];//computing inverse now to multiply, than divide later
        //        }}

        //        {for( int k = 0; k < numk; k++ )
        //        {
        //            kseedsl[k] = sigmal[k]*inv[k];
        //            kseedsa[k] = sigmaa[k]*inv[k];
        //            kseedsb[k] = sigmab[k]*inv[k];
        //            kseedsx[k] = sigmax[k]*inv[k];
        //            kseedsy[k] = sigmay[k]*inv[k];
        //            //------------------------------------
        //            //edgesum[k] *= inv[k];
        //            //------------------------------------
        //        }}
        //    }
        //}

        //===========================================================================
        ///	PerformSupervoxelSLIC
        ///
        ///	Performs k mean segmentation. It is fast because it searches locally, not
        /// over the entire image.
        //===========================================================================
        //void PerformSupervoxelSLIC(
        //    vector<double>&				kseedsl,
        //    vector<double>&				kseedsa,
        //    vector<double>&				kseedsb,
        //    vector<double>&				kseedsx,
        //    vector<double>&				kseedsy,
        //    vector<double>&				kseedsz,
        //        int**&					klabels,
        //        const int&				STEP,
        //    const double&				compactness)
        //{
        //    int sz = m_width*m_height;
        //    const int numk = kseedsl.size();
        //        //int numitr(0);

        //    //----------------
        //    int offset = STEP;
        //        //if(STEP < 8) offset = STEP*1.5;//to prevent a crash due to a very small step size
        //    //----------------

        //    vector<double> clustersize(numk, 0);
        //    vector<double> inv(numk, 0);//to store 1/clustersize[k] values

        //    vector<double> sigmal(numk, 0);
        //    vector<double> sigmaa(numk, 0);
        //    vector<double> sigmab(numk, 0);
        //    vector<double> sigmax(numk, 0);
        //    vector<double> sigmay(numk, 0);
        //    vector<double> sigmaz(numk, 0);

        //    vector< double > initdouble(sz, DBL_MAX);
        //    vector< vector<double> > distvec(m_depth, initdouble);
        //    //vector<double> distvec(sz, DBL_MAX);

        //    double invwt = 1.0/((STEP/compactness)*(STEP/compactness));//compactness = 20.0 is usually good.

        //    int x1, y1, x2, y2, z1, z2;
        //    double l, a, b;
        //    double dist;
        //    double distxyz;
        //    for( int itr = 0; itr < 5; itr++ )
        //    {
        //        distvec.assign(m_depth, initdouble);
        //        for( int n = 0; n < numk; n++ )
        //        {
        //                        y1 = max(0.0,			kseedsy[n]-offset);
        //                        y2 = min((double)m_height,	kseedsy[n]+offset);
        //                        x1 = max(0.0,			kseedsx[n]-offset);
        //                        x2 = min((double)m_width,	kseedsx[n]+offset);
        //                        z1 = max(0.0,			kseedsz[n]-offset);
        //                        z2 = min((double)m_depth,	kseedsz[n]+offset);


        //            for( int z = z1; z < z2; z++ )
        //            {
        //                for( int y = y1; y < y2; y++ )
        //                {
        //                    for( int x = x1; x < x2; x++ )
        //                    {
        //                        int i = y*m_width + x;

        //                        l = m_lvecvec[z][i];
        //                        a = m_avecvec[z][i];
        //                        b = m_bvecvec[z][i];

        //                        dist =			(l - kseedsl[n])*(l - kseedsl[n]) +
        //                                        (a - kseedsa[n])*(a - kseedsa[n]) +
        //                                        (b - kseedsb[n])*(b - kseedsb[n]);

        //                        distxyz =		(x - kseedsx[n])*(x - kseedsx[n]) +
        //                                        (y - kseedsy[n])*(y - kseedsy[n]) +
        //                                        (z - kseedsz[n])*(z - kseedsz[n]);
        //                        //------------------------------------------------------------------------
        //                        dist += distxyz*invwt;
        //                        //------------------------------------------------------------------------
        //                        if( dist < distvec[z][i] )
        //                        {
        //                            distvec[z][i] = dist;
        //                            klabels[z][i]  = n;
        //                        }
        //                    }
        //                }
        //            }
        //        }
        //        //-----------------------------------------------------------------
        //        // Recalculate the centroid and store in the seed values
        //        //-----------------------------------------------------------------
        //        //instead of reassigning memory on each iteration, just reset.

        //        sigmal.assign(numk, 0);
        //        sigmaa.assign(numk, 0);
        //        sigmab.assign(numk, 0);
        //        sigmax.assign(numk, 0);
        //        sigmay.assign(numk, 0);
        //        sigmaz.assign(numk, 0);
        //        clustersize.assign(numk, 0);

        //        for( int d = 0; d < m_depth; d++  )
        //        {
        //            int ind(0);
        //            for( int r = 0; r < m_height; r++ )
        //            {
        //                for( int c = 0; c < m_width; c++ )
        //                {
        //                    sigmal[klabels[d][ind]] += m_lvecvec[d][ind];
        //                    sigmaa[klabels[d][ind]] += m_avecvec[d][ind];
        //                    sigmab[klabels[d][ind]] += m_bvecvec[d][ind];
        //                    sigmax[klabels[d][ind]] += c;
        //                    sigmay[klabels[d][ind]] += r;
        //                    sigmaz[klabels[d][ind]] += d;

        //                    clustersize[klabels[d][ind]] += 1.0;
        //                    ind++;
        //                }
        //            }
        //        }

        //        {for( int k = 0; k < numk; k++ )
        //        {
        //            if( clustersize[k] <= 0 ) clustersize[k] = 1;
        //            inv[k] = 1.0/clustersize[k];//computing inverse now to multiply, than divide later
        //        }}

        //        {for( int k = 0; k < numk; k++ )
        //        {
        //            kseedsl[k] = sigmal[k]*inv[k];
        //            kseedsa[k] = sigmaa[k]*inv[k];
        //            kseedsb[k] = sigmab[k]*inv[k];
        //            kseedsx[k] = sigmax[k]*inv[k];
        //            kseedsy[k] = sigmay[k]*inv[k];
        //            kseedsz[k] = sigmaz[k]*inv[k];
        //        }}
        //    }
        //}


        //===========================================================================
        ///	SaveSuperpixelLabels
        ///
        ///	Save labels in raster scan order.
        //===========================================================================
        //void SaveSuperpixelLabels(
        //    const int*&					labels,
        //    const int&					width,
        //    const int&					height,
        //    const string&				filename,
        //    const string&				path) 
        //{
        ////#ifdef WINDOWS
        //        char fname[256];
        //        char extn[256];
        //        _splitpath(filename.c_str(), NULL, NULL, fname, extn);
        //        string temp = fname;
        //        string finalpath = path + temp + string(".dat");
        //////#else
        //////        string nameandextension = filename;
        //////        size_t pos = filename.find_last_of("/");
        //////        if(pos != string::npos)//if a slash is found, then take the filename with extension
        //////        {
        //////            nameandextension = filename.substr(pos+1);
        //////        }
        //////        string newname = nameandextension.replace(nameandextension.rfind(".")+1, 3, "dat");//find the position of the dot and replace the 3 characters following it.
        //////        string finalpath = path+newname;
        //////#endif

        //        int sz = width*height;
        //    ofstream outfile;
        //    outfile.open(finalpath.c_str(), ios::binary);
        //    for( int i = 0; i < sz; i++ )
        //    {
        //        outfile.write((const char*)&labels[i], sizeof(int));
        //    }
        //    outfile.close();
        //}


        //===========================================================================
        ///	SaveSupervoxelLabels
        ///
        ///	Save labels in raster scan order.
        //===========================================================================
        //void SaveSupervoxelLabels(
        //    const int**&				labels,
        //    const int&					width,
        //    const int&					height,
        //    const int&					depth,
        //    const string&				filename,
        //    const string&				path) 
        //{
        ////#ifdef WINDOWS
        //        char fname[256];
        //        char extn[256];
        //        _splitpath(filename.c_str(), NULL, NULL, fname, extn);
        //        string temp = fname;
        //        string finalpath = path + temp + string(".dat");
        ////#else
        ////        string nameandextension = filename;
        ////        size_t pos = filename.find_last_of("/");
        ////        if(pos != string::npos)//if a slash is found, then take the filename with extension
        ////        {
        ////            nameandextension = filename.substr(pos+1);
        ////        }
        ////        string newname = nameandextension.replace(nameandextension.rfind(".")+1, 3, "dat");//find the position of the dot and replace the 3 characters following it.
        ////        string finalpath = path+newname;
        ////#endif

        //        int sz = width*height;
        //    ofstream outfile;
        //    outfile.open(finalpath.c_str(), ios::binary);
        //    for( int d = 0; d < depth; d++ )
        //    {
        //        for( int i = 0; i < sz; i++ )
        //        {
        //            outfile.write((const char*)&labels[d][i], sizeof(int));
        //        }
        //    }
        //    outfile.close();
        //}

        //===========================================================================
        ///	EnforceLabelConnectivity
        ///
        ///		1. finding an adjacent label for each new component at the start
        ///		2. if a certain component is too small, assigning the previously found
        ///		    adjacent label to this component, and not incrementing the label.
        //===========================================================================
        //void EnforceLabelConnectivity(
        //    const int*					labels,//input labels that need to be corrected to remove stray labels
        //    const int					width,
        //    const int					height,
        //    int*&						nlabels,//new labels
        //    int&						numlabels,//the number of labels changes in the end if segments are removed
        //    const int&					K) //the number of superpixels desired by the user
        //{
        ////	const int dx8[8] = {-1, -1,  0,  1, 1, 1, 0, -1};
        ////	const int dy8[8] = { 0, -1, -1, -1, 0, 1, 1,  1};

        //    const int dx4[4] = {-1,  0,  1,  0};
        //    const int dy4[4] = { 0, -1,  0,  1};

        //    const int sz = width*height;
        //    const int SUPSZ = sz/K;
        //    //nlabels.resize(sz, -1);
        //    for( int i = 0; i < sz; i++ ) nlabels[i] = -1;
        //    int label(0);
        //    int* xvec = new int[sz];
        //    int* yvec = new int[sz];
        //    int oindex(0);
        //    int adjlabel(0);//adjacent label
        //    for( int j = 0; j < height; j++ )
        //    {
        //        for( int k = 0; k < width; k++ )
        //        {
        //            if( 0 > nlabels[oindex] )
        //            {
        //                nlabels[oindex] = label;
        //                //--------------------
        //                // Start a new segment
        //                //--------------------
        //                xvec[0] = k;
        //                yvec[0] = j;
        //                //-------------------------------------------------------
        //                // Quickly find an adjacent label for use later if needed
        //                //-------------------------------------------------------
        //                {for( int n = 0; n < 4; n++ )
        //                {
        //                    int x = xvec[0] + dx4[n];
        //                    int y = yvec[0] + dy4[n];
        //                    if( (x >= 0 && x < width) && (y >= 0 && y < height) )
        //                    {
        //                        int nindex = y*width + x;
        //                        if(nlabels[nindex] >= 0) adjlabel = nlabels[nindex];
        //                    }
        //                }}

        //                int count(1);
        //                for( int c = 0; c < count; c++ )
        //                {
        //                    for( int n = 0; n < 4; n++ )
        //                    {
        //                        int x = xvec[c] + dx4[n];
        //                        int y = yvec[c] + dy4[n];

        //                        if( (x >= 0 && x < width) && (y >= 0 && y < height) )
        //                        {
        //                            int nindex = y*width + x;

        //                            if( 0 > nlabels[nindex] && labels[oindex] == labels[nindex] )
        //                            {
        //                                xvec[count] = x;
        //                                yvec[count] = y;
        //                                nlabels[nindex] = label;
        //                                count++;
        //                            }
        //                        }

        //                    }
        //                }
        //                //-------------------------------------------------------
        //                // If segment size is less then a limit, assign an
        //                // adjacent label found before, and decrement label count.
        //                //-------------------------------------------------------
        //                if(count <= SUPSZ >> 2)
        //                {
        //                    for( int c = 0; c < count; c++ )
        //                    {
        //                        int ind = yvec[c]*width+xvec[c];
        //                        nlabels[ind] = adjlabel;
        //                    }
        //                    label--;
        //                }
        //                label++;
        //            }
        //            oindex++;
        //        }
        //    }
        //    numlabels = label;

        //    if(xvec) delete [] xvec;
        //    if(yvec) delete [] yvec;
        //}


        //===========================================================================
        ///	RelabelStraySupervoxels
        //===========================================================================
        //void EnforceSupervoxelLabelConnectivity(
        //    int**&						labels,//input - previous labels, output - new labels
        //    const int&					width,
        //    const int&					height,
        //    const int&					depth,
        //    int&						numlabels,
        //    const int&					STEP)
        //{
        //    const int dx10[10] = {-1,  0,  1,  0, -1,  1,  1, -1,  0, 0};
        //    const int dy10[10] = { 0, -1,  0,  1, -1, -1,  1,  1,  0, 0};
        //    const int dz10[10] = { 0,  0,  0,  0,  0,  0,  0,  0, -1, 1};

        //    int sz = width*height;
        //    const int SUPSZ = STEP*STEP*STEP;

        //    int adjlabel(0);//adjacent label
        //        int* xvec = new int[SUPSZ*10];//a large enough size
        //        int* yvec = new int[SUPSZ*10];//a large enough size
        //        int* zvec = new int[SUPSZ*10];//a large enough size
        //    //------------------
        //    // memory allocation
        //    //------------------
        //    int** nlabels = new int*[depth];
        //    {for( int d = 0; d < depth; d++ )
        //    {
        //        nlabels[d] = new int[sz];
        //        for( int i = 0; i < sz; i++ ) nlabels[d][i] = -1;
        //    }}
        //    //------------------
        //    // labeling
        //    //------------------
        //    int lab(0);
        //    {for( int d = 0; d < depth; d++ )
        //    {
        //        int i(0);
        //        for( int h = 0; h < height; h++ )
        //        {
        //            for( int w = 0; w < width; w++ )
        //            {
        //                if(nlabels[d][i] < 0)
        //                {
        //                    nlabels[d][i] = lab;
        //                    //-------------------------------------------------------
        //                    // Quickly find an adjacent label for use later if needed
        //                    //-------------------------------------------------------
        //                    {for( int n = 0; n < 10; n++ )
        //                    {
        //                        int x = w + dx10[n];
        //                        int y = h + dy10[n];
        //                        int z = d + dz10[n];
        //                        if( (x >= 0 && x < width) && (y >= 0 && y < height) && (z >= 0 && z < depth) )
        //                        {
        //                            int nindex = y*width + x;
        //                            if(nlabels[z][nindex] >= 0)
        //                            {
        //                                adjlabel = nlabels[z][nindex];
        //                            }
        //                        }
        //                    }}

        //                    xvec[0] = w; yvec[0] = h; zvec[0] = d;
        //                    int count(1);
        //                    for( int c = 0; c < count; c++ )
        //                    {
        //                        for( int n = 0; n < 10; n++ )
        //                        {
        //                            int x = xvec[c] + dx10[n];
        //                            int y = yvec[c] + dy10[n];
        //                            int z = zvec[c] + dz10[n];

        //                            if( (x >= 0 && x < width) && (y >= 0 && y < height) && (z >= 0 && z < depth))
        //                            {
        //                                int nindex = y*width + x;

        //                                if( 0 > nlabels[z][nindex] && labels[d][i] == labels[z][nindex] )
        //                                {
        //                                    xvec[count] = x;
        //                                    yvec[count] = y;
        //                                    zvec[count] = z;
        //                                    nlabels[z][nindex] = lab;
        //                                    count++;
        //                                }
        //                            }

        //                        }
        //                    }
        //                    //-------------------------------------------------------
        //                    // If segment size is less then a limit, assign an
        //                    // adjacent label found before, and decrement label count.
        //                    //-------------------------------------------------------
        //                    if(count <= (SUPSZ >> 2))//this threshold can be changed according to needs
        //                    {
        //                        for( int c = 0; c < count; c++ )
        //                        {
        //                            int ind = yvec[c]*width+xvec[c];
        //                            nlabels[zvec[c]][ind] = adjlabel;
        //                        }
        //                        lab--;
        //                    }
        //                    //--------------------------------------------------------
        //                    lab++;
        //                }
        //                i++;
        //            }
        //        }
        //    }}
        //    //------------------
        //    // mem de-allocation
        //    //------------------
        //    {for( int d = 0; d < depth; d++ )
        //    {
        //        for( int i = 0; i < sz; i++ ) labels[d][i] = nlabels[d][i];
        //    }}
        //    {for( int d = 0; d < depth; d++ )
        //    {
        //        delete [] nlabels[d];
        //    }}
        //    delete [] nlabels;
        //    //------------------
        //    if(xvec) delete [] xvec;
        //    if(yvec) delete [] yvec;
        //    if(zvec) delete [] zvec;
        //    //------------------
        //    numlabels = lab;
        //    //------------------
        //}

        //===========================================================================
        ///	DoSuperpixelSegmentation_ForGivenSuperpixelSize
        ///
        /// The input parameter ubuff conains RGB values in a 32-bit unsigned integers
        /// as follows:
        ///
        /// [1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1]
        ///
        ///        Nothing              R                 G                  B
        ///
        /// The RGB values are accessed from (and packed into) the unsigned integers
        /// using bitwise operators as can be seen in the function DoRGBtoLABConversion().
        ///
        /// compactness value depends on the input pixels values. For instance, if
        /// the input is greyscale with values ranging from 0-100, then a compactness
        /// value of 20.0 would give good results. A greater value will make the
        /// superpixels more compact while a smaller value would make them more uneven.
        ///
        /// The labels can be saved if needed using SaveSuperpixelLabels()
        //===========================================================================
        public void DoSuperpixelSegmentation_ForGivenSuperpixelSize(
            uint* ubuff,
            int width,
            int height,
            int* klabels,
            int numlabels,
            int superpixelsize,
            double compactness)
        {
            ////------------------------------------------------
            //const int STEP = sqrt(double(superpixelsize))+0.5;
            ////------------------------------------------------
            //vector<double> kseedsl(0);
            //vector<double> kseedsa(0);
            //vector<double> kseedsb(0);
            //vector<double> kseedsx(0);
            //vector<double> kseedsy(0);

            ////--------------------------------------------------
            //m_width  = width;
            //m_height = height;
            //int sz = m_width*m_height;
            ////klabels.resize( sz, -1 );
            ////--------------------------------------------------
            //klabels = new int[sz];
            //for( int s = 0; s < sz; s++ ) klabels[s] = -1;
            ////--------------------------------------------------
            //if(1)//LAB, the default option
            //{
            //    DoRGBtoLABConversion(ubuff, m_lvec, m_avec, m_bvec);
            //}
            //else//RGB
            //{
            //    m_lvec = new double[sz]; m_avec = new double[sz]; m_bvec = new double[sz];
            //    for( int i = 0; i < sz; i++ )
            //    {
            //            m_lvec[i] = ubuff[i] >> 16 & 0xff;
            //            m_avec[i] = ubuff[i] >>  8 & 0xff;
            //            m_bvec[i] = ubuff[i]       & 0xff;
            //    }
            //}
            ////--------------------------------------------------
            //bool perturbseeds(false);//perturb seeds is not absolutely necessary, one can set this flag to false
            //vector<double> edgemag(0);
            //if(perturbseeds) DetectLabEdges(m_lvec, m_avec, m_bvec, m_width, m_height, edgemag);
            //GetLABXYSeeds_ForGivenStepSize(kseedsl, kseedsa, kseedsb, kseedsx, kseedsy, STEP, perturbseeds, edgemag);

            //PerformSuperpixelSLIC(kseedsl, kseedsa, kseedsb, kseedsx, kseedsy, klabels, STEP, edgemag,compactness);
            //numlabels = kseedsl.size();

            //int* nlabels = new int[sz];
            //EnforceLabelConnectivity(klabels, m_width, m_height, nlabels, numlabels, double(sz)/double(STEP*STEP));
            //{for(int i = 0; i < sz; i++ ) klabels[i] = nlabels[i];}
            //if(nlabels) delete [] nlabels;
        }

        //===========================================================================
        ///	DoSuperpixelSegmentation_ForGivenNumberOfSuperpixels
        ///
        /// The input parameter ubuff conains RGB values in a 32-bit unsigned integers
        /// as follows:
        ///
        /// [1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1]
        ///
        ///        Nothing              R                 G                  B
        ///
        /// The RGB values are accessed from (and packed into) the unsigned integers
        /// using bitwise operators as can be seen in the function DoRGBtoLABConversion().
        ///
        /// compactness value depends on the input pixels values. For instance, if
        /// the input is greyscale with values ranging from 0-100, then a compactness
        /// value of 20.0 would give good results. A greater value will make the
        /// superpixels more compact while a smaller value would make them more uneven.
        ///
        /// The labels can be saved if needed using SaveSuperpixelLabels()
        //===========================================================================
        public void DoSuperpixelSegmentation_ForGivenNumberOfSuperpixels(
            uint* ubuff,
            int width,
            int height,
            int* klabels,
            int numlabels,
            int K,//required number of superpixels
            double compactness)//weight given to spatial distance
        {
            //const int superpixelsize = 0.5+double(width*height)/double(K);
            //DoSuperpixelSegmentation_ForGivenSuperpixelSize(ubuff,width,height,klabels,numlabels,superpixelsize,compactness);
        }

        //===========================================================================
        ///	DoSupervoxelSegmentation
        ///
        /// There is option to save the labels if needed.
        ///
        /// The input parameter ubuffvec holds all the video frames. It is a
        /// 2-dimensional array. The first dimension is depth and the second dimension
        /// is pixel location in a frame. For example, to access a pixel in the 3rd
        /// frame (i.e. depth index 2), in the 4th row (i.e. height index 3) on the
        /// 37th column (i.e. width index 36), you would write:
        ///
        /// unsigned int the_pixel_i_want = ubuffvec[2][3*width + 36]
        ///
        /// In addition, here is how the RGB values are contained in a 32-bit unsigned
        /// integer:
        ///
        /// [1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1]  [1 1 1 1 1 1 1 1]
        ///
        ///        Nothing              R                 G                  B
        ///
        /// The RGB values are accessed from (and packed into) the unsigned integers
        /// using bitwise operators as can be seen in the function DoRGBtoLABConversion().
        ///
        /// compactness value depends on the input pixels values. For instance, if
        /// the input is greyscale with values ranging from 0-100, then a compactness
        /// value of 20.0 would give good results. A greater value will make the
        /// supervoxels more compact while a smaller value would make them more uneven.
        //===========================================================================
        void DoSupervoxelSegmentation(
            uint** ubuffvec,
            int width,
            int height,
            int depth,
            int** klabels,
            int numlabels,
            int supervoxelsize,
            double compactness)
        {
            ////---------------------------------------------------------
            //const int STEP = 0.5 + Math.Pow(double(supervoxelsize),1.0/3.0);
            ////---------------------------------------------------------
            //vector<double> kseedsl(0);
            //vector<double> kseedsa(0);
            //vector<double> kseedsb(0);
            //vector<double> kseedsx(0);
            //vector<double> kseedsy(0);
            //vector<double> kseedsz(0);

            ////--------------------------------------------------
            //m_width  = width;
            //m_height = height;
            //m_depth  = depth;
            //int sz = m_width*m_height;

            ////--------------------------------------------------
            //    //klabels = new int*[depth];
            //m_lvecvec = new double*[depth];
            //m_avecvec = new double*[depth];
            //m_bvecvec = new double*[depth];
            //for( int d = 0; d < depth; d++ )
            //{
            //            //klabels[d] = new int[sz];
            //    m_lvecvec[d] = new double[sz];
            //    m_avecvec[d] = new double[sz];
            //    m_bvecvec[d] = new double[sz];
            //    for( int s = 0; s < sz; s++ )
            //    {
            //        klabels[d][s] = -1;
            //    }
            //}

            //DoRGBtoLABConversion(ubuffvec, m_lvecvec, m_avecvec, m_bvecvec);

            //GetKValues_LABXYZ(kseedsl, kseedsa, kseedsb, kseedsx, kseedsy, kseedsz, STEP);

            //PerformSupervoxelSLIC(kseedsl, kseedsa, kseedsb, kseedsx, kseedsy, kseedsz, klabels, STEP, compactness);

            //EnforceSupervoxelLabelConnectivity(klabels, width, height, depth, numlabels, STEP);
        }

    }
}
